home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 188_01 / trans.c < prev    next >
Text File  |  1987-09-29  |  7KB  |  217 lines
  1. /*
  2. HEADER:        ;
  3. TITLE:        C elementary transcendentals;
  4. VERSION:    1.0;
  5.  
  6. DESCRIPTION:       Source    code    for   all   standard    C 
  7.                transcendentals
  8.  
  9.         Employs  ldexp()  and  frexp()  functions;  if 
  10.                suitable  versions  of these are  not  provided 
  11.                by  a given compiler,  the versions provided in 
  12.                source  code  will require  adaptation  to  the 
  13.                double float formats of the compiler.
  14.  
  15. KEYWORDS:    Transcendentals, library;
  16. SYSTEM:        CP/M-80, V3.1;
  17. FILENAME:    TRANS.C;
  18. WARNINGS:      frexp()   and  ldexp()  are   implementation 
  19.                dependent.   The  compiler  employed  does  not 
  20.                support    minus   (-)   unary   operators   in 
  21.                initializer  lists,  which are required by  the 
  22.                code.
  23. AUTHORS:    Tim Prince;
  24. COMPILERS:    MIX C, v2.0.1;
  25. */
  26. /* Copyright (C) 1986, Chandler Software,
  27. ** 4456 W. Maple Rd., Birmingham MI 48010-1923
  28. ** (313) 855-4587
  29. ** C adaptation of 1984 FORTRAN code
  30. ** prepared for non-profit distribution by C Users' Group */
  31. /*$NESTCMNT*/
  32. # include "a:stdio"
  33. /* dblfmt describes floating point format
  34. ** machine dependencies in case frexp() and ldexp() are not
  35. ** supplied */
  36. union dblfmt{
  37.   double flt;
  38.   struct{
  39.     char mant[7]; /* full reverse byte order
  40. ** IEEE int sign:1;unsigned expn:11;unsigned mant:4
  41. ** with DEC reverse byte pairing, can't cross 16-bit
  42. ** boundaries */
  43.     char expn;
  44.   }fld;    /* char same as unsigned :8 in this C */
  45. }
  46. #define ROUND(x) (int)((x>=0)-.5+x) /* Pascal ROUND */
  47. #define ODD(i) ((i)&1) /* also like Pascal */
  48. /* teach the preprocessor some algebra */
  49. #define PI 3.1415926535897932385
  50. #define HPI 1.5707963267948966192
  51. #define LN2 .69314718055994530942
  52. #define L2E 1.442695040888963407
  53. #define cos(x) sin(HPI-(x)) /* double sin() */
  54. #define tan(x) sin(x)/cos(x) /* double sin() */
  55. /* for future use
  56. #define fabs(x) (x>=0?x:-(x))
  57. /* double sin() */
  58. #define atan2(x,y) (y>0?atan((x)/(y)):y==0?(x>=0?HPI:-HPI):(x>=0?PI+atan((x)/(y)):atan((x)/(y))-PI))
  59. #define asin(x) atan2(x,sqrt((1-(x))*(1+x))) /* double sin(),sqrt() */
  60. #define acos(x) atan2(sqrt((1-(x))*(1+x)),x) /* double sin(),sqrt() */
  61. #define log(x) log2(x)*LN2 /* double log2() */
  62. #define exp(x) exp2((x)*L2E) /* double exp2() */
  63. #define cosh(x) sqrt(sinh(x)*sinh(x)+1) /* double sinh(),sqrt() */
  64. #define tanh(x) (x<-20?-1.:(x>=20?1.:sinh(x)/cosh(x)))
  65. */
  66. #define pow(x,y) exp2(log2(x)*(y))
  67. main(){ /* test accuracy of tan,sin,cos,atan,exp2,log2 */
  68.   double sin(),atan(),log2(),exp2();
  69.   double x;
  70.   for(x=1e-36;x<1e33;x*=1e4)
  71.     printf("x:%23.16e tan(atan):%23.16e \n\tpow:%23.16e\n",
  72.       x,tan(atan(x)),pow(x,1.));
  73. }
  74. double frexp(x,scale)
  75.   int *scale;
  76.   double x;
  77. {
  78.   union dblfmt x1;
  79.   #define BIAS 128 /* exponent field of .9 (IEEE 1023) */
  80.   x1.flt=x;
  81.   *scale=x1.fld.expn-BIAS; /* scale by .5^(*scale) */
  82.   x1.fld.expn=BIAS;  /* .5 <= result < 1 */
  83.   return(x1.flt);
  84. }
  85. double ldexp(x,scale)
  86.   double x;
  87.   int scale;
  88. {
  89.   union dblfmt x1;
  90.   x1.flt=x;
  91. /* scale by 2^scale */
  92.   x1.fld.expn+=scale;
  93.   return(x1.flt);
  94. }
  95. double sin(x)
  96.   double x;
  97. {
  98.   #define NTS 8
  99. /* negative signs are ignored by some compilers! */
  100.   static double table[NTS]={
  101.     -.7370812e-12,.160478576e-9,-.2505187133e-7,
  102.     .275573164289e-5,-.00019841269823293,
  103.     .008333333333276252,-.1666666666666597,.9999999999999999};
  104.   double poly(),pm,floor();
  105. /* try to take care of x>maxint unless much too slow
  106. ** employ periodicity sin(x+n*PI)=(-1)^n*sin(x) */
  107.   x-=PI*(pm=floor(x/PI+.5));
  108. /* now |x| < HPI; use series and fix sign */
  109.   return(poly(x*x,NTS,table)*(ODD((int)pm)?-x:x));
  110. }
  111. /* use fast polynomial evaluation (possibly non-portable) */
  112. double poly(x,n,table)
  113.   double x;
  114.   int n;
  115.   double table[];
  116. {
  117.   register double sum;
  118.   register int i;
  119.   sum=table[0];
  120. /* unroll loop by pairing terms to save overhead */
  121.   for(i=2;i<n;i+=2)
  122.     sum=(sum*x+table[i-1])*x+table[i];
  123.   return(i==n?sum*x+table[i-1]:sum);
  124. }
  125. double atan(x)
  126.   double x;
  127. {
  128.   #define TNPI6 .57735026918962576451
  129.   #define NTA 9
  130.   static double table[NTA]={
  131.     .0443911883527,-.06483241134718,.07679374240257,
  132.     -.090903714847573,.111110979102203,-.1428571410307564,
  133.     .1999999999873222,-.33333333333329944,
  134.     .99999999999999999};
  135.   char invert,reduce,neg;
  136.   double poly();
  137. /* atan(-x)=-atan(x) */
  138.   if(neg=(x<0))x=-x;
  139. /* tan(HPI-x)=1/tan(x) */
  140.   if(invert=(x>=1))x=1/x;
  141. /* tan(a-b)={tan(a)-tan(b)}/{tan(a)*tan(b)+1}; b=PI/6 */
  142.   if(reduce=(x>=.269))x=(x-TNPI6)/(x*TNPI6+1);
  143.   x*=poly(x*x,NTA,table);
  144.   if(reduce)x+=.52359877559829887308;
  145.   if(invert)x=HPI-x;
  146.   return(neg?-x:x);
  147. }
  148. double log2(x)
  149.   double x;
  150. {
  151.   #define NTL 7
  152.   static double table[NTL]={
  153.     .24291838162,.2614440423316,.3206163946133,.41219840197457,
  154.     .57707801724629,.961796693924339,2.885390081777927};
  155.   int scale,incsc;
  156.   double poly(),ldexp(),frexp();
  157. /* split x -> 2^scale*(reduced x near 1) */
  158.   incsc=frexp(x,&scale)<.7071;
  159.   x=ldexp(x,-(scale-=incsc));
  160.   x=(x-1)/(x+1);
  161.   return(poly(x*x,NTL,table)*x+scale);
  162. }
  163. double exp2(x)
  164.   double x;
  165. {
  166.   #define INFINITY 1.7e38/*IEEE 0x7ff0000000000000*/
  167.   #define MAXEXP 127 /* IEEE 1023 */
  168.   #define MINLG2 -128
  169.   double evenp,oddp,xsq,ldexp();
  170.   int scale;
  171.   if(x>=MAXEXP)return(INFINITY);
  172.   if(x<MINLG2)return(0);
  173. /* 2^x=2^ROUND(x)*2^(x-ROUND(x)) */
  174.   x-=scale=ROUND(x);
  175. /* approximation is ratio of polynomials differing only in
  176. ** sign of odd terms */
  177.   xsq=x*x;
  178. /* continued fraction approx for e^x
  179.   (x*(15120+xsq*(420+xsq))+30240+xsq*(3360+xsq*30))/... */
  180. /* minimax fit for 2^x, 18 significant digits */
  181.   oddp=x*(65556.325877407443+xsq*(874.804294426022+xsq));
  182.   evenp=189155.572484473087+xsq*(10097.515376265486+
  183.     xsq*43.302654542155);
  184.   return(ldexp((evenp+oddp)/(evenp-oddp),scale));
  185. }
  186. /* for future use
  187. double sinh(x)
  188.   double x;
  189. {
  190.   #define NTSH 7
  191.   static double table[NTSH]={
  192.     1.632721e-10,2.50484370e-8,2.7557344083e-6,
  193.     1.984126975507e-4,.0083333333334753,.1666666666666579,
  194.     1.0000000000000001};
  195.   double poly(),exp2();
  196.   return(fabs(x)<1?poly(x*x,NTSH,table)*x:(exp(x)-1/exp(x))/2);
  197. }
  198. double sqrt(x) /* as if division hardware but no sqrt */
  199.   double x;
  200. { static float table[]={.59,.417,.295}; /* minimax 2 digits */
  201.   register double x1;
  202.   double frexp(),ldexp();
  203.   register int i;
  204.   register char iters=3; /* number of Newton iterations */
  205.   register float x0; /* early stages need 2.5 digits precision */
  206.   int scale;
  207.   if(x<=0)return(x);
  208.   x0=frexp(x,&scale); /* sqrt(2^x*y)=2^(x/2)sqrt(y) */
  209.   i=ODD(scale); /* separate fits for y<.5 & y>.5 */
  210.   x1=ldexp(x0*table[i]+table[i+1],(scale+1)>>1); /*crude sqrt*/
  211.   do /* enough Newton iterations for full accuracy */
  212.     x1=ldexp(x/x1+x1,-1); /* ldexp might be faster than *.5 */
  213.   while(--iters);
  214.   return(x1);
  215. }
  216. */
  217.